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ABSTRACT 


An explicit and unconditionally stable finite difference method for the solution 
of the transient inverse heat conduction problem in a semi-infinite or finite 
slab mediums subject to nonlinear radiation boundary conditions is presented. 
After measuring two interior temperature histories, the mollification method is 
used to determine the surface transient heat source if the energy radiation law 
is known. Alternatively, if the active surface is heated by a source at a rate 
proportional to a given function, the nonlinear surface radiation law is then 
recovered as a function of the interface temperature when the problem is 
feasible. Two typical examples corresponding to Newton cooling law and Stefan 
-Boltzmann radiation law respectively are illustrated. In all cases, the 
method predicts the surf ace conditions with an accuracy suitable f or many 
practical purposes. 
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1. Introduction. 

In this paper we investigate the numerical identification of surface transient 
heat sources in one-dimensional semi-infinite and finite slab mediums when the 
active surface radiates energy according to a known nonlinear law. Alternatively 
if the active surface is heated by a source at a rate proportional to a given 
function, the nonlinear radiating boundary condition is then numerically identified 
as a function of the interface temperature if the problem is feasible. 

These two tasks can be viewed as suitable generalizations of the classical 
problem of attempting to determine the interface temperature between a gas and a 
solid with a nonlinear heat transfer law. The existence and uniqueness of a 
strictly increasing solution of the semi-infinite body version of this problem has 
been considered by Mann and Wolf [Ref.7] for a monotone Lipschitz radiation law. 
Roberts and Mann [Ref. 10] extended the previous result after removing the Lipschitz 
condition on the nonlinear heat transfer law. Keller and Olmstead [Ref. 6] 
investigated the same problem in the presence of a positive integrable transient 
source and introduced a constructive proof for existence and uniqueness of the 
interface temperature by the method of lower and upper solutions. The numerical 
solution of the nonlinear Volterra integral equation characterizing the active 
surface temperature history was implemented by Chambre [Ref.l] using the method of 
succesive approximations and, more recently, by Groetsch [Ref. 3] who succesfully 
combined Abel inversion formula with B-spline approximation and product 
integration. A natural extension of this technique to solve the same problem in the 
finite slab medium is discussed in Groetsch [Ref.4], Also for the finite slab case, 
Villasenor and Squire [Ref. 12] have proposed a numerical procedure based on a 
generalized trapezoidal rule and Richardson extrapolation. More general problems of 
the same kind, combining the effects of convection and radiation at the interface, 
can be found in Friedman [Ref. 2] and Saljnikov and Petrovic [Ref. 11]. 

In all the works mentioned above, the nonlinear radiation law and the transient 
boundary source are supposed to be known in order to determine the interface 
temperature. Consequently, if the new task consists on the identification of the 
nonlinear radiation law or on the identification of the transient boundary source 
function, a different approach must be used. 

It is possible to estimate the surface temperature and the surface heat flux in 
a body from measured temperature histories at fixed locations inside the body 
However, this Inverse Heat Conduction Problem (IHCP) is an ill-posed problem 
because small errors in the data induce large errors in the computed surface heat 
flux history or in the computed temperature history solutions and, consequently, 
special methods are needed in order to restore continuity with respect to the data. 
In this paper we consider initially, the solution of a one-dimensional IHCP by a 
fully explicit and stable space marching finite difference implementation of the 
Mollification Method introduced by Murio [Ref. 8] and Guo, Murio and Roth [Ref. 5]. 
The procedure allows for a direct discretization of the differential equation and 
it is generated by automatically filtering the noisy data by discrete mollification 
against a suitable averaging kernel and then using finite differences, marching in 
space, to numerically solve the associated well-posed problem. Once the temperature 
and the heat flux transient functions have been approximately recovered at the 
interface, it is a simple task to numerically identify the transient heat source if 
the nonlinear radiation law is known. On the other hand, if the surface is heated 
by a source at a rate proportional to a given function, we proceed to approximately 
recover the nonlinear surface radiation law describing the physical conditions at 
the interface, provided that the range of temperatures at the interface contain 
sufficient information. 

In Section 2, we define the new identification problems with data specified on a 
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p continuum of time and data errors measured in the L 2 norm and derive rigorous 
stability bounds. The efficiency of the method is demonstrated in Section 3, where 
together with a description of the numerical procedure, we present the results of 
several computational experiments with rapidly varying and discontinuous profiles, 
for both linear - Newton cooling law - and nonlinear - Stefan-Boltzmann law - 
models. In all cases, numerical stability and good accuracy are achieved even for 
small time steps and high levels of noise in the data. Section 4 includes a summary 
and some conclusions. 

2. Description of the Problem. 

We consider a one-dimensional IHCP in a semi-infinite or finite slab, in which 
the temperature and heat flux histories f C t ) and q(t) on the left-hand surface (x = 
0) are desired and unknown, and the temperature and heat flux at some interior 
point x = xq or at the right-hand surface x = a are approximately measurable. Note 

that, equivalently, the data temperature histories might be measured at two 
interior points. For the semi-infinite medium, 0 < xq and for the finite slab, 0 < 
xq - a. We assume linear heat conduction with constant coefficients and normalize 
the problem by dimensionless quantities. Without loss of generality, we consider xq 
= a = 1 in all cases. The problem can be described mathematically as follows. 

For the semi-infinite or finite slab, the unknown temperature u(x,t) satisfies 
respectively, 


u t (x,t) = U xx (x,t), 

t > 0, 

u(l.t) = F(t), 

t > 0, 

-U x (l,t) = Q(t), 

t > 0, 

u(x,0) = u 0 (x), 

0 < x 

u(0,t) = f(t). 

t > 0, 

-u x (0,t) = q(t) 



= E(u(0,t))-g(t), t > 0, 


0 < x < oo or 0 < x < 1, 

with corresponding approximate 
data function Fm^> 
with corresponding approximate 
data function 
<oo or 0 < x < 1, 

the desired but unknown 
temperature function, 

the desired but unknown 
heat flux function. 


(la) 

(lb) 

(lc) 

(l d) 

(le) 


(If) 


The nonlinear boundary condition (If), indicates that the active surface radiates 
energy at a rate proportional to E and is heated at a rate proportional to the 
function g. Our aim is to obtain more detailed information about the boundary 
condition at the interface x = 0. More precisely, we want to estimate the 
function E if g is known or, reciprocally, we want to identify the source function 
g if the radiation law E is given. 

We also assume that all the functions involved are L 2 functions in any time 
interval of interest and use the corresponding L 2 norm, as defined below, to 
measure errors: 

Ilf II = [ s\ Z | f (t) | 2 dt ] 1/2 . 

In this setting, it is also natural to hypothesize that the exact data functions 
F(t) and Q(t) and the measured data functions F m (t) and Q m (t) satisfy the L 2 data 
error bounds 


L 


II F— F m If ^ c 


and 


HQ-Q m ll * c. 


j 
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f lt is weI1 ^own that solving for f(t) and q(t) from F(t) and Q(t) amplifies] 

every Fourier frequency component of the error by the factor exp[w/2] 1/2 , -oo <w< oo. 
This shows that the inverse problem is highly ill-posed in the high frequency 
components. See Murio [Ref. 8] and Guo, Murio and Roth [Ref. 5] for further 
discussions. 

Stabilized Problem . 

The one-dimensional IHCP can be stabilized if instead of attempting to find the 
point values of the temperature function f(t) or the heat flux function q(t), we 
attempt to reconstruct the 6-mollification of the functions f and q at time t, 
given by 

J 6 f(t) = (P 6 *f)(t), J 5 q(t) = (p 6 «q)(t), 

where 

pJt) = exp[-t 2 /6 2 ] 

6 6 * 1/2 

is the one-dimensional Gaussian kernel of radius 5 > 0. The mollifier p (t) is 

5 

always positive, falls to nearly zero outside the interval centered at the origin 
and radius 36 and 


00 

(p^*f)(t) = J p^(x)f(t-T) dr 

-oo 


is the one-dimensional convolution of the functions p and f. We notice that J f(t) 

o 6 

is a C 00 (infinitely differentiable) function and that the mollifier has total 
integral 1. Mollifying system (1), we obtain the following associated problem: 
Attempt to find = J^u(O.t) and J£q m (t) = -J^u x (0,t) at some point t of 

interest and for some radius 6 > 0, given that J u(x,t) satisfies for the semi- 

o 

infinite or finite slab respectively, 


( J 6 u) t ^ J 5 u ^xx» 
J 5 u(1,t) = 

J^ulx.O) = J g u 0 (x,0), 
J 6 u(°,t) = J 6 f m (t), 
-J s u x (0,t) = J 6 q m (t), 


t>0, 0<x<ooor0<x<l, 
t > 0, 
t > 0, 

0<x<coor0<x<l, 
t > 0, unknown, 
t > 0, unknown. 


( 2 ) 


This problem and its solutions satisfy the following: 

Theorem 1. Suppose that IIF-F m ll ^ c and IIQ-Q m il ^ c. Then 

(i) Problem (2) is a formally stable problem with respect to perturbations in the 
data. 

(ii) If the exact boundary temperature function f(t) and the exact heat flux 
function q(t) have uniformly bounded first order derivatives on the bounded domain 
D = [0,T], then J^f m and J^q m verify 
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[“ Hf-JgfmUD - °( 5 > + 3c exp[(26)' 2/3 ] (3) | 

and 

llq-J^qm^D - + - (1 + 3 exp[6" 2/3 ]). (4) 

The proof of this statement can be found in Guo, Murio and Roth [Ref. 5]. 

Once the mollified temperature and mollified heat flux functions have been 
evaluated at the interface, it is feasible to attempt to identify the source 
function g or the radiation energy function E given in formula (If). 

Identification of the source f unction g. 

Assuming that the radiation law at the active surface is known, according to 
(If), the exact source function is given by 

g(t) = E(f(t))-q(t). (5) 

The approximate source function, denoted g a (t), is defined by 

g a (t) = E(J 6 f m (t))-J g q m (t), (6) 

and in order to estimate the error, we suppose that the surface radiates energy at 

a rate proportional to [f(t)]*\ Here p is a positive integer, the value p = 1 
corresponding to Newton's law of cooling and p = 4 to Stefan’s radiation law. 

The difference (5) - (6) gives 

g(t) - g a (t) = [f(t)] P - [J g f m (t)] P + q(t) - J g q m (t). 

From the identity a n - b n = (a-b)(a n_ Ca n- ^b +...+ab n-2 +b n- M, taking norms and 
introducing M = max {II J^fmUco.D’ llpll oo,D^> we g et 

llg - g a ll D s P M P_1 Ilf - J g f m ll D + ii q - 

Combining the last inequality with the upper bounds (3) and (4), we obtain the 
estimate 

Hg - ga n D s (pM P_1 + 1) {0(6) + 3c exp[5- 2/3 ]}. (7) 

This shows that the identification of the source function g is stable with respect 
to errors in the data functions F and Q , for fixed p and 5 > 0. 

Remarks : 

1. Notice that the approximate source function g a is actually a function of the 

radius of mollification 5, the amount of noise in the data c and the exponent p in 
the radiation model E. 

2. From a more theoretical point of view, inequality (7) can be used to show the 
convergence of g a to g in the L 2 norm. In fact, setting 0(6) = C 6 for some 

constant C > 0, and choosing 6 = [ln(l/e 1/2 )]” 3/2 , after replacing these quantities 
in (7), we obtain 

I'g - ga^D 5 (pM P_1 + l)<C[ln(l/e 1/2 )]- 3/2 + 3c ,/2 >. 

This last inequality implies that, for the special selection of the radius of 
mollification indicated above, llg 11 D ^ 0 as e 0, for any value of p. 

Identification of the radiation law f unction E. 




L 


From equation (If) it follows that the exact function E, assuming that the 
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T source function g is given, satisfies 

1 E(u(0,t)) = E(f (t)) = g(t) + q(t). (8 ) 

The approximate function, denoted E a , is defined by 

Ea^/rntt)) = + (9) 

Subtracting (8) from (9), taking norms and using inequality (4), we inmediately 
have 


HE “ Ea^D 5 (1+ 3 exp[6 _2/3 ]). 


This estimate also shows that the identification of the radiation law - as a 
function of time - is stable with respect to perturbations in the data functions F 
and Q, for a fixed 8 > 0, provided that the source function is known. However, this 
information is clearly not sufficient to identify the physical process at the 
interface. Nevertheless, since at each time tj we know the ordered pairs 

(tpJ^fm^i)) a nd (tj,E a (tj)), it is possible to collect the coordinates 

^g^m^i^^a^i^ f° r t in a discrete subset of D and obtain a graph of the 

approximate funcional relationship between the radiation law and the temperature at 
the interface. This is certainly always the case if the cardinality of the range of 
temperatures is sufficiently large. Similar remarks to the ones in the 

previous paragraph, about the parameter dependency of E a and convergence in the L 2 
norm of E a to E as the quality of the data functions improve, e ^ 0, also apply 
here. 

The computational details are presented in the next section. 

3. Numerical Procedure. 


With v - J^u and z - -dv/3x, system (2) is equivalent to 
dv dz 

~ ~ > t>0, 0<x<ooor0<x<l, 

ot dx 


dv 

2 “ t>0, 0<x<ooor0<x<l, 

v(l,t) = J g F m (t), t > 0, 

z(l,t) = J 5 Q m (t), t > 0, (11) 

v(x,0) = J 5 u 0 (x,0), 0 < x < 00 or 0 < x < 1, 

v(0,t) = Jgf m (t), t > 0, unknown, 

z(0,t) = JgQ m (t), t > 0, unknown. 

Without loss of generality, we will seek to reconstruct the unknown mollified 
boundary temperature function J g f m and the mollified boundary heat flux function 


J g q m in the unit interval I = [0,1] of the time axis (x = 0). Consider a uniform 


grid in the (x,t) space:{(xj = ih, t n = nk), i = 0,1 N, Nh = 1; n = 0,1 M, 

Mk = L), where L depends on h and k in a way to be specified later, L > 1. 

Let the grid functions V and W be defined by 

v ? = v f x i*t n ), w" = zfxj.tn), OiisN, Osn^M. 

Notice that 


L 
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f vft = J 5 F m (t n ), w£ = J a Q m (t n ), 0 ^ n ^ M, 


and 

v ? = J^otxj.O), 0 i i s N. 

We approximate the partial differential equation in system (11) with the consistent 
finite difference schemes 


w?., = w? - | ( vr‘ 

v ?-l = V? - h W{_ lf 



( 12 ) 


i = N, N-l,...,l; n = 1,2,...,M-1. 


Notice that, as we march backward in the x-direction, we must drop the estimation 
of the interior temperature from the highest previous point in time. Since we want 

to evaluate {Vq} and {Wq} at the grid points of the unit time interval I = [0,1] 

after N iterations, the minimum initial length L of the data sample interval in the 
time axis (x = 1) needs to satisfy the condition L = kM = 1 - k +k/h. 

Once the temperature J^f m an ^ the ^ ux h ave been reconstructed, we 

proceed with the approximate identification of the source function g a or the 
radiation law function E a as explained in Section 2. 


Remarks : 

1. The radius of mollification, 5, can be selected automatically as a function of 
the level of noise in the data. In fact, for a given c > 0, there is a unique 6 > 
O, such that 

N'JgFm - F m ll D = e. (13) 

For the proof of this assertion and some discussions on the numerical 
implementation of this practical selection criterion, see Murio [Ref. 9]. 

2. For the proof of the unconditional stability of the finite difference scheme 
(12) and the analysis of the convergence of the numerical solution of the mollified 
problem (11), the reader should consult Guo, Murio and Roth [Ref. 5]. 

Numerical Results. 


In order to test the accuracy and the stability properties of our method, in 
Problem 1, the approximate reconstruction of a source function g(t) and a nonlinear 
radiation law E(u(0,t)) are investigated for a one-dimensional finite slab exposed 
to a heat flux data function at the free surface x = 1 given by -u x (l,t) = Q(t) = 

0, t > 0, and a temperature data function 


L 


u(l,t) = F(t) 


1 2 oo (-l)n 

(t-0.2)- - - — £ exp [-n^7r^(t-0.2)l, t > 0.2, 

6 n 2 n =i n 2 

0, 0 < t ^ 0.2. 


The exact source solution to be approximately reconstructed at the interface x = 0 


has equation g(t) = E(u(0,t)) - q(t), where E(u(0,t)) = [u(0,t)] P and q(t) = 

-u x (0,t). We consider the values p = 1 and p = 4 corresponding to Newton’s law of 
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I” cooling and Stefan s radiation law respectively. The exact radiation law at the"! 

I interface is given by E(u(0,t)) = g(t) + q(t) and we only consider the nonlinear ' 

case p = 4. If the initial temperature distribution u(x,0) is zero, the exact 
temperature and heat flux functions at the interface are given respectively by 
a function at the free surface x = 1 given by -u x (l,t) = Q(t) = 0, t > 0, and a 
temperature data function 


r 


u(0,t) = f(t) 


and 


1 2 oo 1 

• ( t-0.2)+ y . — 

2 Tt 2 n=1 n 2 

l °- 


exp [-n 2 7i 2 (t-0.2)] > t > 0.2, 
0 < t £ 0.2. 


-u x (0,t) 



t > 0.2, 
0 < t £ 0.2. 


With this information we generate the exact functions E(u(0,t)) and git) for our 
model problem. 

In Problem 2, we attempt to approximately reconstruct the transient source 
function git) for a semi-infinite body initially at zero temperature with data 
functions 


uil.t) = Fit) 


and 


erfc[(t-0.2)" 1/z /2], t > 0.2, 

0, 0 < t £ 0.2, 


[ n(t-0.2))' 1/2 exp{-[4)t-0.2)l' 1 ), 

0 . 

The unique temperature solution at the interface is 

1, t > 0.2, 

0, 0 < t £ 0.2. 

and the corresponding heat flux at the interface is 


-u x (l,t) = Q(t) = J 


u(0,t) = fit) = 


t > 0.2, 
0 < t £ 0.2. 


— u x (0,t) — q(t) 


j [ 7r(t-0.2)]~ 1/2 , t > 0.2, 

0, 0 < t £ 0.2. 


In this case, we do not attempt the identification of the radiation law at the 
active boundary. The energy as a function of the interface temperature is either 0 
or 1 for any value of p making its identification impossible. There is no enough 
information in the range of boundary temperatures which in this example is reduced 
to just two temperature values. 

Since in practice only a discrete set of points is generally available, we shall 
assume that the data functions F m and Q m are discrete functions measured at equally 
spaced points in the time domain I = [0,L], where L = 1 - k + k/h, Nh = 1 h = Ax 
and k = At. In order to compute J 5 F m (t n ) and J^Q m (t n ) in I, we need to extend the 

data functions in such a way that F m and Q m decay smoothly to zero in the interval 
*Sm&x ” ^~36 majc ,L+36 max ] and both are zero in R - I .In what follows, we 

°max 

consider the extended discrete data functions F m and Q m defined at equally spaced 
sample points on any interval of interest in the time axis. 

The selection of the radius of mollification is implemented by solving the 


J 


L 
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discrete version of equation (13), using the bisection method. 

Once the radii of mollification 5p and <5q, associated with the data functions F m 


and Q m respectively, and the discrete filtered data functions J 5 F m (t n )=VN and 

JgQm^n^WN’ O-n-M, are determined with 6 = maxfSp.Sg), we apply the 

finite difference algorithm described previously in this section, marching backward 

in the x-direction. The values Vq and Wq, 0 s n s M-N, so obtained, are then taken 

as the accepted approximations for the interface temperature and heat flux 
histories respectively at the different time locations at x = 0. Finally, we 
identify the approximate transient source function g a or the approximate radiation 

law function E a at the grid points of the time interval I = [0,1] using equations 
(6) and (9). 

In all cases, we use h = Ax = 0.01 and k = At = 0.01. Thus, N = 100, L = 1.99, M 
= 200, 6 max = 0.1 and I = [-0.3,2.29]. The noisy data is obtained by adding a 

^max 

random error to the exact data at every grid point t n in L : 

^ o 

max 

^m^n^ = F(t n ) + c n j 
Qm^n^ = + c n,2> 

where e n j and e n> 2 are Gaussian variables of variance c 


If the discretized computed transient source function component is denoted by g n 

3 . 

and the true component is g n = g(t n ), we use the sample root mean square norm to 

measure the error in the discretized interval I = [0,1], The solution error is then 
given by 


»g a ^ gi'l 


1 M-N 

- — - M I (g a 

M-N n = 1 a 


g n > 2 



If the discretized computed radiation law function component is denoted by E a = 
E a (t n ) and the true component is E n = E(t n ), after evaluating the ordered pairs 

0 - n - M-N, we obtain a graph of the approximate functional relationship 

between the radiation law and the temperature at the interface. This plot is then 
compared with the exact graph corresponding to the values (f (t n ),E(t n )) of the 
model problem. 

Tables 1 and 2 show the results of our numerical experiments associated with 
Problems 1 and 2 respectively, when attempting to identify the transient source 
function at the interface. In all cases, the numerical stability of the method is 
confirmed. The uniformly smaller error norms in Problem 1 are expected since 
at time t = 0.2 the exact source solution has a finite jump discontinuity while in 
Problem 2 the exact source solution has an infinite jump at time t = 0.2. For this 
reason, we have added an extra column in Table 2 indicating the error norms in the 
time interval [.3,1], after the discontinuity. It is clear that the method rapidly 
dissipates the effect of the singularity, a very desirable feature. 

The qualitative behavior of the reconstructed transient source function for 
Problem 1 is illustrated in Figures 1 and 2 where the numerical solution for an 
average perturbation c = 0.005 (full line) is plotted for p = 1 (Newton’s cooling 
law) and p = 4 (Stefan’s radiation law) respectively. In Figure 3 we show the graph 
associated with the reconstructed nonlinear radiation law as a function of the 
approximate temperature at the interface for p = 4 (full line) and the exact 

L J 
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P boundary radiation law (star symbols). Figures 4 and 5 show the computed source 
functions (full lines) for p = 1 and p = 4 respectively, for Problem 2 and for the 
noise level c = 0.005. 


PROBLEM 1 

P 

= 1 (Newton ) 

p = 4 ( Stef an ) 

c 

5 

Error norm 

e 

6 

Error norm 

0.000 

0.04 

0.0866 

0.000 

0.04 

0.0867 

0.002 

0.06 

0.0921 

0.002 

0.06 

0.0929 

0.005 

0.06 

0 . 1014 

0.005 

0.06 

0.1038 


Table 1. Error norm as a function of the level of noise 


j| PROBLEM 2 ~j| 

P = 

= 1 ( Newon ) 

p = 4 (Stefan) 

e 

6 

Error norm 
10, 1 ]/[ . 3 , 1 ] 

G 

6 

Error norm 
[0,1 )/[ .3, 1 ] 

0.000 

0.04 

.5208/. 0183 

0.000 

"0.04 

.5135/. 0373 

0.002 

0.06 

.556 0/. 063 1 

0.002 

0.06 

.5673/. 0675 

0.005 

0.06 

.5879/. 1108 

0.005 

0.06 

.5935/. 1375 


Table 2. Error norm as a function of the level of noise 


4. Conclusions. 

An explicit and unconditionally stable space marching finite difference method 
for the solution of the one-dimensional transient inverse heat conduction problem 

has been implemented for the numerical identification of surface heat sources, if 

the energy radiation law at the active interface is known, and to the numerical 

identification of the nonlinear surface radiation law if the surface is heated by a 
source at a rate proportional to a given function and the interface temperature 
contains enough information. 

The computational procedure is applied to two examples corresponding to Newton 
cooling law and to Stefan— Boltzmann radiation law. In both problems, the source 
functions to be identified have discontinuous histories and in one case an infinite 
jump. The algorithm restores stability with respect to the data, which is essential 
for the introduction of the inverse problem approach, and good accuracy is 

obtained, even for small time sample intervals and relative high noise levels in 
the data. 
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Fig. 2 Source function for Stefan Law 
Problem 1 , 0=0.005, 8=0.06, At=0.01 

Exact: (* * *) ; Computed: ( ) 
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Fig. 3 Reconstructed Stefan radiation Law 
Problem 1 , 0=0.005, 8=0.06, At=0.01 

Exact: (* * *) ; Computed: ( ) 
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Fig. 4 Source function for Newton Law 
Problem 2, 0=0.005, 8=0.06, At=0.01 

Exact: (* * *) ; Computed: ( ) 


Fig. 5 Source function for Stefan Law 
Problem 2, 0=0.005, 8=0.06, At=0.01 

Exact: (* * *) ; Computed: ( ) 



